## This Rfile creates "Figure 2. Number of FNA operators (bars) and FNA income by status (lines)"
##  File written by Kazuma Takakura kazuma.arukakat@gmail.com (University of Tokyo)
##  Last updated on April 27, 2024


setwd("/Users/kazuma.takakura/Dropbox/Saguling_FN_study/data") ## change directory where the data locates
library(ggplot2)
library(dplyr)
library(tidyr)
second_rate = 0.1

df_plot <- read.csv("figuredata.csv") %>%
  mutate(fnnumscaled = fnnum*second_rate) 

# 各軸の範囲をきめる
# yも最小値〜最大値の形で設定してよいが明示的に与えたい場合が多い気がする
x.lim  <- c(min(df_plot$yr), max(df_plot$yr))
y1.lim <- c(0, 140)
y2.lim <- c(0, 7)

# 左軸と右軸の関係を表すスケーラをつくる
# 各軸の最大最小差の比をとっている
scaler <- (y1.lim[2] - y1.lim[1])/(y2.lim[2] - y2.lim[1])

df_plot %>% 
  ggplot() + 
  geom_bar(aes(x=yr, y=fnnum), stat='identity', fill="white", colour = "black",width=0.7) +
  geom_line(aes(x=yr,y=fnincome_status0*scaler, linetype = "Exit before 2000"), size=1) +
  geom_line(aes(x=yr,y=fnincome_status1*scaler, linetype = "Continue after 2000"), size=1) +
  geom_line(aes(x=yr,y=realincomemean*scaler, linetype = "Whole"), size=1) +
  scale_y_continuous(limit=y1.lim, expand = c(0, 0), 
                     sec.axis=sec_axis(trans = ~ ./scaler, 
                                       name="Income from FNA (million)")) +
  theme(legend.position="bottom")+
  labs(x="Year", y="Number of active FNA operators", linetype = "FNA Status", color = "")+
  theme_bw()+
  theme(legend.position="bottom")
ggsave("figure_income.png")
